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Abstract 

O ' Combinig the harmonic balance method (HBM) and a continuation method is a 

well-known technique to follow the periodic solutions of dynamical systems when a 
^ , control parameter is varied. However, since deriving the algebraic system containing 

I the Fourier coefficients can be a highly cumbersome procedure, the classical HBM 

Ph ■ is often limited to polynomial (quadratic and cubic) nonlinearities and/or a few 

. harmonics. Several variations on the classical HBM, such as the incremental HBM or 

the alternating frequency/time domain HBM, have been presented in the literature 
C/^ ' to overcome this shortcoming. Here, we present an alternative approach that can 

Q . be applied to a very large class of dynamical systems (autonomous or forced) with 

(-H \ smooth equations. The main idea is to systematically recast the dynamical system in 

quadratic polynomial form before applying the HBM. Once the equations have been 
rendered quadratic, it becomes obvious to derive the algebraic system and solve it by 
the so-called ANM continuation technique. Several classical examples are presented 
to illustrate the use of this numerical approach. 
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1 Introduction 



In the field of engineering, especially when dealing with nonlinear vibrations, 
it is often required to compute the periodic solutions of a system of nonlinear 
differential equations ( [T|2] ). In the following, we focus on numerical meth- 
ods to find periodic solutions. They are usually subdivided into two main 
groups : those relying on the time-domain formulation and those relying on 
the frequency-domain formulation. The first group consists of methods where 
a time integration algorithm, which is generally limited to a single period, is 
used to transform the original differential system into a system of algebraic 
equations, which are then solved by continuation. With these methods, the 
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unknowns in the algebraic system are the values of the original unknown vari- 
ables at grid points along the periodic orbit. The classical shooting technique 
( [3f4] ) and the orthogonal collocation methods used in AUTO ([5ii7|) belong 
to this group. The second group corresponds to the so-called harmonic bal- 
ance method where the unknown variables are decomposed into truncated 
Fourier series. In this case, the unknowns in the final algebraic system, which 
is obtained by balancing the harmonics featuring in the differential equations, 
are simply the Fourier coefficients of the original unknowns. Methods of both 
kinds are widely used in many applications and research is still required to 
improve them. In practical terms, the choice between the time-domain or the 
frequency-domain approach depends mainly on whether the periodic solution 
can be decomposed with a few Fourier components. The method chosen to 
deal with a given problem can therefore depend on the operating point. 

In this paper, we focus on the harmonic balance method (HBM) and some 
of its variations, which are briefly recalled here. Although the name "har- 
monic balance" seems to have flrst appeared in 1936 ([8]), this method has 
been widely used only since the sixties, and especially for electrical and me- 
chanical engineering purposes. Forced vibrations were flrst studied ([9j), and 
self-sustained oscillations a decade later ((IDI)- The study by Nakhla & Vlach 
([TT]) is often said to be a milestone in the modern formulation of HBM. The 
classical HBM is simple in its principle, but it can be cumbersome or even im- 
practicable, as stated by Peng et al. ([I2]) when the system contains complex 
nonlinearities and when a large number of harmonics is required, ie, more than 
5 or 10. The flrst problem which arises is how to derive the algebraic system 
for the Fourier coefflcients and the second problem is how to solve this strongly 
nonlinear system efflciently. To overcome these shortcomings, many variations 
on the basic HBM have been proposed in the literature. Here, we will simply 
give an overall picture of incremental harmonic balance (IHB) and alternat- 
ing frequency/time domain harmonic balance (AFT) method. With the IHB 
method ([I3]), the incremental-iterative method used for the continuation is 
closely combined with the harmonic principle, so as to be able to apply the HB 
principle to the incremental linear problem instead of the nonlinear original 
system. Nonlinearity of all kinds can be treated using this technique [T4fT5] . 
In the AFT variant, the harmonic balance of the Fourier coefflcients is not 
explicitly performed. At each increment (or iteration) in the continuation pro- 
cedure, the unknowns are transferred to the time domain (using the inverse 
Fast Fourier Transform, which is denoted IFFT below), so as be able to use 
the original system of equations. The nonlinear responses are then transferred 
back to the frequency domain using the Fast Fourier Transform, (denoted 
FFT below). The use of FFT/IFFT procedures seems to have started in |T6] . 
This technique has been used, for instance, to analyse nonlinear vibrations in 
mechanical systems with contact and dry friction [17]. Since the pioneering 
work on IHB and AFT, many variations have been proposed to extend the 
applications, as well as to decrease the computational cost (see fT8] or fT9] for 
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some recent examples). 

An alternative strategy is proposed here for applying the classical HBM with 
a large number of harmonics, without avoiding explicit balancing. The idea 
is to apply a transformation to the original system of differential equations 
before applying the HBM. The aim of this procedure is to transform the non- 
linearities present in the original system into purely polynomial quadratic 
terms. The classical HBM procedure can then be easily applied. This trans- 
formation might seem to be a limitation of the method, since not every system 
can be recast in polynomial quadratic form. However, it will be shown in this 
paper that a very large class of systems with smooth equations can indeed 
be recast in quadratic form by making a few algebraic manipulations and 
a a few additions of equations and auxiliary variables. The idea of using a 
quadratic formulation to simplify the (Fourier) series expansion was inspired 
by another numerical method, the so-called Asymptotic Numerical Method 
(ANM) ( [2T1I22] ! which was presented for the first time in 1990 ([20]). ANM is 
a continuation technique ([23]) involving high order powers series expansion 
of the branches of solutions. For reasons that will be given later on, this ANM 
continuation is the ideal means of solving the algebraic system resulting from 
the HB method proposed here. 

The paper is organised as follows: details of the transformation into a polyno- 
mial quadratic form, the application of the HBM method and the principle of 
continuation are given in section 2 in the cases of autonomous systems. Three 
simple examples are given as an illustration. Section 3 is devoted to periodi- 
cally forced systems, and two further examples are presented. Section 4 gives 
the results of all the examples and ends with some concluding comments and 
perspectives. 



2 Presentation of the method for obtaining periodic solutions of 
autonomous systems 

Let us consider an autonomous system of differential equations 

Y = f{Y,X) (1) 

where 1" is a vector of unknowns, / a smooth nonlinear vector valued function 
and A a real parameter. The dot stands for the derivative with respect to time 
t. We assume that this system has branches of periodic solutions when A varies, 
and we want to find and follow them by applying the harmonic balance method 
and a continuation procedure (path following technique) . The important case 
of a forced (non-autonomous) system will be dealt with separately in section 

El 
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2.1 The harmonic balance principle 

As recalled in the Introduction, the HBM consists basically in decomposing 
Y{t) into a truncated Fourier series : 



This ansatz is put into Eq. ([T]) and f(Y, A) is expanded into Fourier series. By 
balancing the first 2H + 1 harmonic terms, one obtains an algebraic system 
with 2H +1 vector equations for the 2H +1 vector unknowns 1^, the unknown 
pulsation u, and the parameter A. Adding a phase condition ([3]) yields a 
system with equations and A^ + 1 variables. The branches of solutions of 
this algebraic system are then followed using a continuation technique. This 
procedure provides only approximate periodic solutions, since in the expan- 
sion of Eq. all the harmonic terms greater than H remain unbalanced. 
However, if the number of harmonics H is large enough, accurate solutions 
can be obtained. 

The most crucial point is the expansion of the vector f(Y, A) into Fourier 
series. This can be quite a straighforward procedure if f{Y, A) is a polyno- 
mial of degree two or three, and if the number of harmonics H is sufficiently 
small. In this case, the expansion can be carried out by hand. But in the most 
general situation, where f{Y, A) shows nonlinearities of any kinds, this com- 
putation can be very cumbersome, even with the help of a symbolic software 
program. The only alternative is then to use numerical procedures to estimate 
the Fourier series of f(Y, A), by performing iterative FFT/IFFT steps, for ex- 
ample. This drawback was the starting-point in the developpment of many 
variants of the HBM as mentioned in the Introduction. 



2.2 A key point: the quadratic recast 

The main aim of this paper is to present a simple but powerful procedure that 
overcomes the drawback mentioned above. The idea is to recast the original 
system ([T]) into a new system where the nonlinearities are at most quadratic 
polynomials . The application of the harmonic balance method subsequently 
becomes quite straighforward. This new quadratic system, which will be writ- 
ten as follows 



H H 



Y{t) = l^o + E cos{kut) + Ys,k sm{kut) 



(2) 



k=l k=l 



m{Z) =c + l{Z) + q{Z, Z) 



(3) 
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can contain both differential and algebraic equations. A^e is taken to denote the 
number of equations. The unknown vector Z (size A^e) contains the original 
components of the vector Y and some new variables which are added to get 
the quadratic form. The right hand side of Eq. (I3|) is written as follows: c 
is a constant vector with respect to the unknown Z, /(.) is a linear vector 
valued operator with respect to the vector entry, and .) is a quadratic 
vector valued operator, which is linear with respect to both entries. At this 
stage, the three vectors c, /(.) and .) may depend on the real parameter 
A. However, it is preferable for the quadratic operator .) not to depend on 
A, which may require introducting a new variable in Z (see example 1 below), 
on the left hand side, m(.) is a linear vector valued operator with respect to 
the vector entry. The algebraic equations correspond to zero values of m(.). 

Since the nonlinearities are quadratic, the HBM can obviously be easily and 
systematically applied on Eq. ((31), even with a large number of harmonics. The 
question is now: can any system be easily put into the form of Eq. ((31) ? As 
we will see, the answer is yes, for a very large class of nonlinear systems. 

The idea of recasting a nonlinear system into a new one with a quadratic poly- 
nomial nonlinearity is frequently used when considering the so-called ANM 
continuation method, which consists in computing power series expansion of 
solution branches, see f [2T]l22ll23]l24| ) for the basic algorithms and f [25]l26ll27] ) 
for some applications in mechanics. The quadratic recast gives a very simple 
and systematic algebra for the series determination, even with high orders of 
truncature. Similar benefits are obtained with the Fourier series expansions. 

It is now worth illustrating this main idea by giving some elementary exam- 
ples where the recast from Eq. (P) to Eq. (JS]) will be performed explicitely. 
We will deal first with the classical Van der Pol oscillator, the Rossler sys- 
tem and a model for clarinet-like musical instruments. Other examples, with 
nonlinearities of others kinds such as rational fraction, will be given later on. 



Example 1: The Van der Pol oscillator. We take the following second order 
autonomous system, known as the Van der Pol oscillator, where the parameter 
A governs the amplitude of the nonlinear damping term. 

ii- A(l -M^) M + M = (4) 



This equation can be classically recast into first order system ( as Eq. ([T])) by 
introducing the velocity vif) = u{t) as a new unknown. We obtain 



u = V 

V = A(l — u^) V — u 



(5) 
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If we now introduce the auxiliary variables w{t) = 1 — u^(t) and r(t) = 
v{t)w{t), the system can be written in the form of Eq. ((31) with Z = [u, v, w, r]*. 
The various terms are arranged below so that it will be clear to the readers 
how the different functions m, c, I and q are formed. 

ii = +v 

V = —u + Ar 

0=1 -w -u" 

^-2-' ~ +r —v w 

m{Z) ^ 1{Z,X} q{Z,Z) 



We finally obtain two first order differential equations and two algebraic equa- 
tions with quadratic polynomial nonlinearities. Only the operator / depends 
here on A. 



Example 2: The Rossler system. We take the following first order autonomous 
system of three equations known as the Rossler system. 

X = —y — z 

y = x + ay (7) 
z = b + z{x — A) 



where x, y, z are functions of time and a, b, c are float parameters. This 
system can be written in the form of Eq. (JS]) without any auxiliary variables 
(i.e. Z = [x,y,zY): 




(8) 



Here again, only the operator / depends on A. This example will be used 
here to show that the approach presented in this paper can be used to find 
period-doubling and bifurcations with period 4. 

Example 3: A model for clarinet-like musical instruments. In the case of 
small amplitude oscillations, a simple model for reed instruments (clarinet. 
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saxophone, etc) can be written using a modal formulation (see [28]): 

X + qrUJrX + LO^X = UJ^p 

pn + 2anCPn + OO^Pn = f ^ Vn G [1 . . . A^] 

u = c(i — A + x)^yx — p 

P = En=lPn 



where the unknowns x, Pn=i..N, P, u are functions of time and Qr, cOr, a„, a;„, c, 
/, C are given parameters describing either the physics or the player's action, 
is the number of acoustic modes. A stands here for the blowing pressure. 

It is worth noting that even if this system contains a square root, it can be 

recast in the form of Eq. ([3]): if we introduce the auxiliary variables y = x, Zn = 

Pn and V = y/X — p, system (l9|) can be rewritten, with Z = [x, y,pi, . . . ,Pn, zi . . . z^, u, f]* 

as follows: 



X = 





+y 


y 





+tOj.p — QrCOrV — iO^X 


Pn = 





~\~Zn 


Zyi 





-2anCZn - UjIPu 








-u + C{l-X)v 








-p + Pi + ...+Pn 


= ^ 


_A_ 


-p 


m(Z) 




l(Z,X) 



Vn e [1 . . . AT] 



-Qxv 



q{Z,Z) 



Here, the linear term /, but also the constant term c, are functions of A. 



2.3 The harmonic balance method applied to a quadratic system 



In this section, the harmonic balance method is applied to the system (l3|). 
The unknown (column) vector Z is decomposed into Fourier series with H 
harmonics: 

H H 

Z{t) = Zo + ^ Zc^k cos^kut) + ^ Zs^k sm^kujt) (11) 

k=l k=l 
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The components of the Fourier series are collected into a large (column) vector 
[/, with size {2H + 1) x A^e, where Ne is the number of equations in ((31). 




(12) 



By introducing the expansion (fTTl) into the set of Eqs. (I3|), collecting the terms 
of the same harmonic index, and neglecting the higher order harmonics, one 
obtains a large system of {2H + 1) x A/'e equations for the unknown vector f/. 



The new operators M(.), C, L{.), and Q(., •) that apply to U depend only 
on the operators m(.), c, /(.) and g(., .) of Eq. (I3|) and on the number of 
harmonics H. The explicit formulas have been given in annex 1 for a sake of 
brevity. It should be noted that these expressions for M(.), L(.) and Q{.,-) 
are the same in all the examples presented. The change from one particular 
system to another only involves the operators m(.), /(.) and g(., .). 

The final system (fT3l) contains {2H + 1) x A^e equations for the {2H + 1) x A^e 
unknowns U plus the angular frequency u and the continuation parameter A. 
Since the original system is autonomous, a phase condition has to be added 
to Eq. f fT3l) to define a unique orbit. Indeed, the time t does not appear in Eq. 
dS]) and if U{t) is a solution of Eq. (fT3|l then U{t + r) is also a solution, for 
any r. We refer here to textbooks dealing with periodic solution continuation 
( [3l|7l|4| ) for the choice of the phase condition. 

2.4. The continuation procedure 

2.4-1 Framework 

As the result of the harmonic balance operation, we now have to solve an 
algebraic system 



with R e and U = [U\X,uY e M^+^ (A^ = (2i/ + 1) x A^^ + !)• Any 
numerical continuation method ( [29f30f3] ) and any suitable software (see ( [3T] 
for an overview) can be used for this purpose. However, in this paper we will 
use the so-called Asymptotic-Numerical-Method (ANM), on which the idea 
of the quadratic recasting was based. It is worth noting that the final system 
(fT3ll is already quadratic with respect to U and w, and that the application of 
the ANM is therefore quite straighforward. 



uM{U) = C + L{U) + Q(U, U) 



(13) 



R(U) = 



(14) 
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In the continuation procedure, a pseudo-arc length parametrization will be 
used in order to be able to pass the limit points with respect to A. The pa- 
rameter A will therefore now become an unknown, like U and u. It is now 
mandatory to specify how the operators m, c, / and q in Eq. (I3l) depend on A. 
In the three examples presented here, we have organised the terms of equations 
so that the parmeter A only appears in the operators c and /. In addition, we 
have managed to make this dependence linear. Once again, this formulation 
is not limited to the three example presented here, but it can be obtained for 
a very large class of dynamical systems provided suitable additional variables 
are introduced. We recall that the parameter A should not appear in the op- 
erator q, for which the HBM algebra is the more complex. For example, if 
there is a term Xu"^ in Eq. ([3]), it should be rewritten as Xv with the auxiliary 
variable v = u^, and put into I instead of q. 

In what follows, it is therefore assumed that c and / can be written: 

c =.„ + Ac, ^^^j 



where cq, ci, Iq and h are independent of A. Under this assumption, the oper- 
ators C and L simply become 

L(.) =M-) + ALi(0 



Note that the expression for Co, Ci and Lq, Li in terms of cq, ci and Iq, h are 
exactely the same as those given in annexe 1 for C and L in terms of c and /. 
The final algebraic system (fT4l) becomes 

R(U) = LO + L(U) + Q(U,U) (17) 

with U = [U\X,uY and 
LO = Co 

L(U) = LoiU) + ACi (18) 
Q(U,U) = Q{U,U) + XLi{U) - ujM{U) 

In Eq. (fTTll . LO is a constant vector, L(.) is a linear vector valued operator 
and Q(., .) is a bilinear vector valued operator. We will now briefly review the 
ANM continuation technique for solving quadratic system EqfflTl). 
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2.4-2 The ANM continuation 

One of the main particularities of the ANM is that it gives access to branches of 
solution in the form of power series. Assuming that we know a regular solution 
point Uo, the branch of solution passing through this point is computed in 
the form of a power series expansion (truncated at order n) of the pseudo- 
arclength path parameter a = (U — Uo)*Ui, where Ui is the tangent vector 
at Uq: 

U(a) = Uo + aUi + a'Us + a^Ug + . . . + a"U„. (19) 

The series (fTOll is replaced in Eq. (fTTl) and each power of a is equated to zero, 
giving a series of linear systems : 

• order : LO + L(Uo) + Q(Uo,Uo) = 0, which is obvious since Uq is a 
solution of Eq. (flTl) . 

• order 1 : L(Ui) + Q(Uo, Ui) + Q(Ui, Uq) = 0, which can also be written 
Ju()Ui = where Ju„ G is the jacobian matrix of R evaluated at 
Uo. 

• order 2 < p < n : JuoUp + Q(U,, Up_,) = 

The original nonlinear problem has therefore been reduced to a series of n 
linear systems of Ne equations. However, at each order, the linear systems are 
under-dimensioned since they have + 1 unknowns. The additional equations 
required are obtained by inserting Eq. ( flQl ) into the definition of the path 
parameter a given above. This gives: 

• order 1 : U^^Ui = 1 

• order 2 < p < n : U^^Up = 

2.4-3 Comments: 

• The original nonlinear system of N^. equations has been transformed into n 
linear systems of equations, which have to be solved successively (as in 
classical perturbation methods), i.e.. Up is deduced from terms occuring at 
lower orders. 

• The range of utility of a truncated power series is generally limited because 
the series have a finite radius of convergence. Once each Up (p G [1 . . . A^]) 
has been found, the range of utility of the series expansion (fT9l) is defined 
by the value Umax such that 

VaG [0 a™,,], II R(U(a))| I <e, (20) 

where is a user-defined tolerance parameter (see ( [24f32] ) for details of 
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the ttmax calculation). Note that the range of utility Umax is generally ap- 
proximately equal to the radius of convergence of the series (|33j). 
The series expansion and the associated range of utility a^ax only define 
one part of the branch of solutions, which is called a section. To determine 
the entire branch, it is necessary to restart the whole series calculation from 
successively updated starting points Uq. The simplest way of performing 
this updating is to take the end point of a previously computed section as 
the starting point for the next section. This is the principle of continuation 
with the ANM. Finally, the entire branch is given as a succession of different 
branch sections, where the length of each section is automatically given by 
its range of utility amax- 

The ANM continuation approach has the following advantages: 

■ the solution branch is known analytically, section by section. 

■ since all the linear systems to be solved have the same Jacobian matrix, 
the computational cost of the series ( fT9l) is low. 

• since the size of the section is given by the convergence properties of the 
current step, ie, by the values of amax, the ANM continuation algorithm 
does not require any special step-length control strategies. 

■ the ANM continuation algorithm is highly robust, even when the branch 
contains sharp turns. Branch switching can therefore be easily performed, 
using small pertubations in the system. 



2.5 Implementation in the MANLAB software program 



2.5.1 A brief presentation of MANLAB 

MANLAB is an interactive software program for the continuation and bifur- 
cation analysis of algebraic systems, based on ANM continuation. The latest 
version is programmed in Matlab using an objet-oriented approach ([34]). 
MANLAB has a Graphical User Interface (GUI) with buttons, on-line inputs 
and graphical windows for generating, displaying and analysing the bifurca- 
tion diagram and the solution of the system. To enter the system of equations, 
the user has to provide three vector valued Matlab functions corresponding 
to the constant, linear and quadratic operators LO, L and Q. As an example, 
take the biochemical reaction system used by Doedel et Al in (^) 



ri(Mi, «2, A) = 2iii -U2 + 100^:^:^ 
r2(wi, M2, A) = 2ii2 - til + 100^^^ - (A + /i) = 0. 



A =0 

(21) 



Introducing the following additional variables Vi = Ui + V2 = U2 + U2t = 
and ^4 = Y+I^5 the system can be rewritten with quadratic polynomial 
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nonlinearities as follows, 






+ 


2Ui — U2 — 


A 


+ 


lOOwiVs 


= 


-fl 


+ 


2U2 — Ul — 


A 


+ 


100^2^4 


= 





+ 


Vi - Ul 






ul 


= 





+ 


V2 - U2 






ul 


= 


-1 


+ 


V3 




+ 


VIV3 


= 


-1 


+ 


V4 




+ 


V2V4 


= 


LO 




L(U) 






Q(U,U) 





(22) 



The first two equations correspond directly to Eq ([2T|). The last four define 
the auxilliary variables Vi,V2,Vs,V4. The constant, linear and quadratic terms 
in each equation have been split to define the three operators LO, L and Q. In 
the MANLAB software program, the unknowns are collected in a single vector 
U = [ui,U2,Vi,V2,V3,V4, X], and the three operators are given by the following 
vector valued functions ( = 0.05 in the case of this example ) : 



function [LO] 



LO 



function [L] = L(U) 



function [Q] = Q(U,V) 



LO=zeros (6,1) ; 



05: 



L=zeros (6,1) ; 

L(1)=2*U(1)- U(2) -U(7); 

L(2)=2*U(2)- U(l) -U(7); 

L(3)=U(3)-U(1); 

L(4)=U(4)-U(2); 

L(5)=U(5); 

L(6)=U(6); 



Q=zeros (6,1) ; 
C|(1)=100*U(1)*V(5) 
Q(2)=100*U(2)*V(6) 
Q(3)= -U(1)*V(1) 
0(4)= -U(2)*V(2) 
Q(5)= U(3)*V(5) 
5(6)= U(4)*V(6) 



2.5.2 Implementation of the periodic solution continuation in MANLAB 



To compute the branches of periodic solutions, the system has to be first 
recasted in the form of Eq. (JSj) with account of the additional decomposition 
( ITSll . This is probably the most unusual and difficult task for a new user. 
Subsequently, it is only necessary to provide the MANLAB software program 
with the operator m(.), Cq, Ci, /(.) and (?(.,.) . The functions LO, L and 
Q, which are the actual input for MANLAB, have been programmed once for 
all, using the expression in Eq. (fTSll and the formulas given in annex 1. The 
examples presented in this paper are available online ([Gj). 



12 



3 The case of a periodically forced system 



We now focus on periodically forced (non-autonomous) systems: 

Y = f{t,Y,X) (23) 

where / is periodic in t, with the period T (forcing period). We look for periodic 
solutions (responses) with a period pT or -, where p is an integer. A classical 
strategy consists in transforming Eq. (l23l) into an augmented autonomous 
system ([3]), by adding an oscillator with the desired forcing period in the 
system of equations, for instance ([7j). 

In what follows, a direct approach will be used. It consists in expanding the 
forcing term into harmonics, and taking them into account in the balance 
of individual harmonics. This appraoch is illustrated below with two further 
examples. 

Example 4: Forced Duffing oscillator. 

The normalized forced Duffing oscillator is the non-autonomous equation: 

M + + M + = /cos(At) (24) 



We take the damping coefficient /x and the force amplitude / constant, and use 
the forcing angular frequency as the varying parameter A. By using v{t) = u{t) 
and w{t) = v?{t), this equation can be recast as follows 



u = V 

V = /cos(At) —2iiv — u —uw 

= w —uu 

m{Z) c(t,A) l(Z) g(Z,Z) 



(25) 



where Z = [u,v,wY, and the forcing term is deliberately put into c. The 
forcing frequency is now related to the response frequency by putting X = uj 
or possibly, X = puj {p is an integer). The term c(t) is then expanded into 
harmonics with respect to u. 

This results in slight changes in the procedure: 

• because of the synchronization of the response and the forcing, the phase 
condition has to be removed. Note that the parameter A is no longer an 
unknown, since it was chosen as a multiple of uj. In comparison with the case 
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of an autonomous system, both the number of equations and the number of 
unknowns have decreased by one. 
• in the final system (fTTl) (fTSl l . the operators Ci and Li disappear and the 
forcing amphtude / has to be accounted for in Lq 

Lastly, we take the Raylegh-Plesset equation, which is used to model the large 
amplitude vibrations of a gas bubble in a fluid. The forcing is handled slightly 
differently and this example also shows how to cope with a power of —3. 

Example 5: The forced Rayleigh-Plesset equation ([35]). 

Let R be the radius of the vibrating bubble, and Rq the radius at rest. The 
equation of motion of the bubble is 

RR+lte = A{{?^f -l} + Bcos{Xt) (26) 
/ R 

where A and B are flxed. We introduce the normalized radius u = S- and the 

Ho 

(normalized) velocity v = Dividing Eq (|26|) by R^ and defining a = 
b = we get the first order system 

ii = V 

(27) 

uv = a{u ^ — 1) — + bcos{Xt) 



We now introduce the following auxiliary variables x = ^, y = x'^ , and z = v"^, 
and arrive at 

ii = V 

V = a{y'^ — x) — \xz + hx cos{Xt) 



After undergoing this transformation, the system is quadratic but the forcing 
term has now been multiplied by the unknown function x, and it is no longer 
a constant term with respect to the unknown. We introduce another auxilliary 
variable r{t) = cos(At), and replace the term bxcos{Xt) by bxr. Finally, we 
obtain the following system, where the first two equations stand for Eq ( l28l ) 

14 



(28) 



and the last four define the auxilliary variables. 



ii = 





+ 


V 


+ 


i) = 





+ 


{—ax) 


+ 


= 


1 


+ 





+ 


= 





+ 


y 


+ 


= 





+ 


z 


+ 


^ = 


— cos(At) 




r ^ 


+ 


m{Z) 


c{t) 




1{Z) 





ay"^ — ^xz + bxr 
(—ux) 



[—xx) 
[-vv) 



(29) 



q{Z,Z) 



Here the unknown vector is Z = [u, v, x, y, z, r]*, and the forcing term is clearly 
present in the operator c. 



4 Numerical results on selected examples 



In the following selected examples, the numerical results obtained using the 
approach presented in this paper are either compared with time-domain sim- 
ulations to show the validity of our approach, or used to illustrate particular 
features: the influence of the number of harmonics in the case of the Van der 
Pol oscillator, the ability to follow period-doubling bifurcations in that of the 
Rossler system, the ability to follow a direct or inverse Hopf bifurcation in 
that of the clarinet model, and to illustrate a forced system, in that of the 
Duffing oscillator. 



One period of the unknown u in the time domain 



One period of the unknown u in the time domain 



time (s) 

Attractor in the phase space (u,v) 



time (s) 

Attractor in the phase space (u,v) 



























- 










- 




















1 Solution Irom time domain ODE-solver 



Figure 1. The Van der Pol oscillator (Eq. Q): Comparison between the harmonic 
balance results and those obtained with a time domain simulation. Influence of the 
number of harmonics H taken into account ((a) H=10, (b) H=50). 
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Example 1: The Van der Pol oscillator. Numerical results obtained on the 
Van der Pol oscillator (Eq. (Jlj)) are presented in Fig. [H Good agreement was 
obtained with the results of a time domain simulation, performed using Matlab 
ODE solvers. As expected, the number of harmonics H in the solution sought 
by the harmonic balance method had to be adapted to span the bandwith of 
the solution obtained by direct integration. In the case of Fig. [T], A = 3, which 
requires having a relatively high H to match the reference solution. 

Example 2: The Rossler system. The ability to follow period-doubling bi- 
furcations is illustrated in Fig. O Following these bifurcations is more difficult 
in the case of autonomous systems than in forced systems, since the period 
of oscillation is also unknown. Moreover, since a T-periodic solution is also a 
2T-periodic solution, one cannot expect the value of the period obtained by 
the harmonic balance process to be doubled when crossing a period-doubling 
bifurcation. The strategy used here therefore consists in introducing complex 
subharmonic amplitudes as additional unknowns. To detect K period-doubling 
bifurcations, the solution is then sought as: 

H H ^ 

Z{t) = Zo + Y^ Zc^k cos{-j^ujt) + Zs,k sm{—ujt) (30) 

fc=l k=l 

When the solution belongs to the T-periodic solution branch, ^c,fc|fc=i k-i ~ ^ 
and K-i ~ 0- When the solution belongs to the 2r-periodic solution 

branch, Zc^k\k=i..K-i,k=^2K-i = and Zs,k\k=i..K-i,k^2K-i = 0- One practical 
consequence in terms of the computational cost is that in order to span the 
same bandwidth, the number of harmonics has be mulitplied by 2^. This 
case is illustrated in Fig. [2] where if = 10 from the T-periodic solution to 
if = 2^ X 10 = 40 after two period-doubling bifurcations. 

Example 3: The clarinet model. 

Two typical examples of the bifurcation diagrams obtained with the clarinet 
model are shown in Fig. [3l These pictures focus on the oscillation threshold 
in order to show the robustness of the method, even at singular points. As 
can be seen from both pictures, the radius of convergence of the power series 
expansion decreases when approaching a bifurcation point (and hence the 
length of each section of a branch lying between two consecutive points in 
Fig. [3] decreases) . This behaviour is typical of ANM ([33j). 

Example 4: The forced Duffing oscillator 

Figure 4 shows the frequency-amplitude diagram of the response obtained with 
the method presented here. A classical bent resonance curve can be observed, 
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(a) 

Configuration space (x,y,2) 2D projection in the (x,y) plane 
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(c) 



Configuration space (x,y,z) 




2D projection in the (x,y) plane 




- Harmonic balance 
Time domain simulation 
2D projection in the (x,z) plane 2D projection in the (y,z) plane 



X y 

Figure 2. The Rossler system (Eq. jT]) with a = b = 0.2): Comparison between the 
results obtained using the harmonic balance method and a time domain simulation. 
These three figures show the ability of the method to follow the sub-harmonic cascade 
from (a) A = 2.5 (perdiod T), to (b) A = 3.5 (period 2T after a first period-doubling 
bifurcation) and (c) A = 4 (period 4T afbSr a second period-doubling bifurcation) 



y 



Y 



Figure 3. The clarinet model (Eq. JH) with qr = 0.01, ujr = 6597, H = 3, N = 5, 
c = 340, C = 0.2, It = 5.6e - 08, = 4e - 08, Cpy = 1.4, r = 0.007, wi = 2426.5, 
W2 = 7279.5, W3 = 12132.5, = 16985.4, = 21838.4, ai = 0.11, 02 = 0.19, 
as = 0.25, = 0.3, as = 0.34): depending on the length of the resonator, the 
bifurcation will be either a direct Hopf bifurcation (left picture, L = 0.22) or an 
inverse Hopf bifurcation (right picture, L = 0.2). 




0.5 1 1.5 2 2.5 3 0.1 0.2 0.3 0.4 0.5 



Frequency, CO Frequency, to 

Figure 4. Duffing oscillator, Eq. ([Ml with /x = 0.1 and / = 1.25 as in (p]). The 
figure shows the amplitude of the harmonics 1, 3 and 5 versus the frequency forcing 
LJ. Solid line : results when 5 harmonics are included (H=5 in Eq. (fTTI) ). Dash line 
: results when 9 harmonics are included. The figure on the right is a zoom of the 
region uj < 0.6. 



as well as some additional peaks corresponding to superharmonic resonances. 
Note that only the individual amplitude Ai = {ul^ + u^^) of the odd har- 
monics have been plotted, since the even harmonics {Aq, A2, A^, . . .) are zero. 
The computation of this branch of periodic orbits required 25 steps of MAN- 
continuation when 5 harmonics were included, and 35 steps when 9 harmonics 
were included. The curves obtained for harmonics 1 and 3 were slightly hanged 
when shifting from H = 5 to H = 9 harmonics in Eq. ( fTTll . The curve obtained 
for A5 is more affected, which confirms the logical result that more than 5 har- 
monics have to be used to obtain an accurate result with harmonic 5. We also 
checked that the curve obtained for A^ was only very slightly modified when 
shifting from H = 9 to H = 11 harmonics in Eq. ( fTTl l. 
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5 Discussion and Conclusion 



The key idea in this study is the quadratic recast. It is worth pointing out 
that the present method contains two quadratic recasts. The first one, which 
occurs in the time domain, transforms the original system ^ into ((31), and 
is presented in section 12. 2[ The second one transforms the algebraic system 
( fTSll obtained with the HBM, into the final system ( flTll which is solved by the 
ANM. The operators in ^ and ^ are linked by ((161) and ([I8]). This point 
accounts for many of the characteristics of the method presented here: 

• The whole procedure takes place in the frequency domain using an analyt- 
ical expression for the system to be solved (as with the classical harmonic 
balance) and allows to use an arbitrarily high number of harmonics (as with 
the AFT procedures). 

• In the ANM continuation, no iterative process are used to solve the algebraic 
system, and hence no problems of convergence arise. The computational cost 
is almost predictible and fairly rather low. On the contrary, convergence of 
the iterative process can be difficult with AFT procedures ([39l sect. 4.6, 
pll7]). 

• There is no need for FFT/IFFT procedures, which increases the computa- 
tional load in the case of AFT approaches. 

• No temporal discretization has to be performed (since the whole procedure 
is in the frequency domain), and there is therefore no need to cope with 
aliasing. 

• There is no need to use numerical schemes to compute derivatives. Since the 
problem is a quadratic one , the Jacobian is calculated analytically and this 
can be done very easily. This is a considerable advantage of the approach 
presented here. Indeed, as stated by [39l pl7], usually "analytical calculation 
of the Jacobian involves considerable calculation and transformation" but 
on the other hand, when considering numerical calculation of the Jacobian, 
"for very nonlinear circuits, the error introduced in the Jacobian estima- 
tion results in inaccurate updates of the unknowns and a large number of 
iterations and possibly no convergence". 

As stated in |T2], to apply HB, "it is always necessary to write specific com- 
putation programs for different nonlinear models". Here, on the contrary, the 
procedure is simple and rather general because of the natural splitting between 
the automated part, which is independent of the problem (i.e. the transforma- 
tion from operators c,l,q to operators LO, L, Q and the ensuing solving using 
the ANM), and the part which depends on the particular problem to be solved 
(the writing of operators c,l,q) which is let to the users. From the point of 
view of the users, this means that the work required is relatively easy and the 
software is simple to use. In addition, very few parameters have to be defined 
by the users: the number of harmonics in the solution and parameter in Eq. 
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(l20ll . which controls the accury of the ANM computation. 

In this paper, we have presented only small-sized systems for the sake of sim- 
plicity. The method can be obviously extended in principle to large systems of 
equations although this may require some additional technical work, and in- 
troduce limitations due to the lack of computer ressources. In [36], the method 
has been applied to compute the forced responses of geometrically nonlinear 
elastic structures discretized using the finite element method. It is well known 
that with problems of this kind, the equations of motion are second order 
differential equations in terms of the displacement nodal variables u, and that 
they show a polynomial cubic nonlinarity. Choosing the nodal velocity v and 
the second Piola-Kirchhoff S (at the Gauss point) as auxiliary variables, the 
system of equations can easily be written in the form of a first order sys- 
tem with quadratic polynomial non linearities. When a cosine forcing term is 
applied, they read (in the discrete form) 

it = V 

Mv = - Jy{BL + BNL{u)ySdv + Fcos{Xt) (31) 
S = D{Bl + \Bnl{u))u 

where M is the mass matrix, D is the Hooke operator and Bl and Bj^l are the 
linear and nonlinear strain-displacement operators |37|. In this form, the gov- 
erning equations consist in a group of first order differential equations (velocity 
definition, equation of motion) and a group of algebraic equations (constitu- 
tive law at the Gauss point of the finite elements). Both groups are quadratic 
in u and S. The algebra given in annex 1 has been directely introduced into 
a home-made FEM code, at the element level. 

To conclude, it is proposed to discuss the limitations of this method and 
suggest some possible lines of future research. First of all, it is worth noting 
that when the original system is not directly in the same quadratic form as 
Eq. ((31), introducting auxiliary variables, as done in this paper, examples may 
lead to additional (possibly non-physical) solutions. For example, in the case 
of the clarinet model (example 3), a variable v = — p is introduced, and 
then further defined by the quadratic equation f ^ = 7 — p. This definition 
corresponds to f = ±a/7 — p. In practical terms, this means that the user has 
to check that the branch of solution computed corresponds to a positive v. 

The main limitation of the method presented is that the original systems 
cannot all be put into the quadratic form Eq. (I3|). For this reason, problems 
where one or several relations are given in the frequency domain (typically, 
when some variables are linked through an impedance) cannot be treated. 
Similar difficulties arise when the original system contains functions (trigono- 
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metric functions, exponential functions, non-integer powers function, etc) that 
are applied to the unknowns. The classical pendulum model 6 + sm{9) = is 
a good example of systems of this kind. In these cases, we can generally intro- 
duce new variables and new differential equations, so that the functions are 
generated by the system itself. In the case of pendulum model, v = 9 and the 
two new variables s{t) = sin(6'(t)) and c(t) = cos(6'(t)) are introduced. Taking 
the derivative with respect to time for the last two equations, we obtain 



V 

s 
c 



V 

—s 

cv 
—sv 



(32) 



This quadratic system corresponds to the pendulum model, provided the fol- 
lowing two initial conditions are added 

s(0) = sin(^(0)) ^^^^ 
c(0) = cos(6'(0)) 



The application of the HBM to Eq. (l32ll will now be quite straighforward, 
but, because of Eq. (l33l) . the final algebraic system will not be quadratic. In 
these cases, the application of the ANM continuation is still feasible but this 
requires more elaborate algebra for computating the power series, as explained 
in ( |38|32| ). The efficiency of the present method on the pendulum model will 
be tested in a near future. Another interesting idea would be to test whether 
this method can be used to solve regularized non-smooth dynamic problems, 
and for instance, the case of vibrating systems with contact conditions and 
friction laws. 
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A Appendix 



Here, we give the expression for the operators M, C, L and Q in Eq. (fT3ll in 
terms of m, c, /, q and H. We recall that Eq. (fT3l) is obtained by substituting 
Eq. (fTTi l into Eq. ((31) and by collecting the terms with the same harmonic 
index (cosines and sines). The inputs of M, L and Q are the vector U, which 
contains the Fourier coefficients of Z(t) 



A.l Constant term 

The constant vector c is associated with the harmonic zero. The operator C 
is simply 



A. 2 Linear term 

For the linear operator /, we have 

l{Z{t)) = l{Zo + J2k=i ^c,k cos{kujt) + Zs^k sm{kujt)) 

= 1{Zq) + /(^c,i) cos(wt) + /(^s,i) sin(a;t) + l{Zc^2) cos(2wt) + . . . 

The operator L is 



L{U) = [1{Z^)\ l{Zs,i)\ 1{Z,,2Y, l{Zs,2)\ 1{Zc,h)\ 1{Zs,hY]' (A.4) 




(A.l) 



C = [c*,0*,0*,...]* 



(A.2) 



For the linear operator m, we have 



m{Z{t)) 



m(J2k=i Zs,kkuJ cos{kujt) — Zc^kkuj sin(A;a;t)) 
+ m{Zs^i)uj cos{ujt) — m{Zc^i)uJ sin{ujt) 
+2m{Zs^2)^ cos(2u;t) + . . . 



(A.5) 



The operator M is 



M{U) = [0*,m(Z,,i)*,-m(Z,,i)*,2m(Z,,2)*,-2m(Z,,2)*,... 
...,Hm{Z,^Hy,-Hm{Z,,HyY 



(A.6) 
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A. 3 Quadratic term 



With the notation Z^^^ — Zq and — Zc^iCos{ujt) + jsin(a;t), the decom- 
position of Z(t) is conveniently rewritten as 



H 



(A.7) 



i=0 



For the operator q, we have 

q{X{t),Y{t)) = Ef=o Ef=o 

= g(X(o),F(o)) +^^^^^(^W^y(o)) (A.8) 



When i > 1 and j > 1, q{X^^\Y^^^) give harmonics i + j and z — j as follows: 

g(XW,yO)) = i{5(X,,,,y,,,.) - q{X,,,,Y,j)}cos{{i+j)u;t) 

mq{X,^i,Ys,j) + q{Xs,^,Yc,J)}sm{{i + j)ujt) 
+ i{g(Xc,j, Ycj) + q{Xs,u Ys,j)] cos((i - j)ixit) 
+ \{(l{Xs,i, yc,j) - q{X^,i, y^j)} sin((i - j)ujt) 



By grouping the terms with the same harmonics index, and canceling any 
harmonics higher than index H, we have 

H 

q{X{t),Y{t)) = go + X] ^c,it cos(A;a;i) + sin(/i;a;t) (A. 10) 

k=l 

with 

go = g(Xo, yo) + E \{q{X,,, Y,,) + q{Xs,, Y,,)} (A.ll) 
and when i > 1 

H 

^1 1 (A.12) 



23 



J2 li-l(^c,j,Y,, 



(A.13) 



j=i+l 
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